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1. Worm algorithms 

Monte-Carlo simulation algorithms based on all-order strong coupling or hopping parameter 
expansions have been known for a long time Yet only the relatively recent ideas of Prokof 'ev 
and Svistunov [Q] opened the way to the formulation of highly efficient algorithms. Instead of 
sampling the partition function of a model the new class of algorithms samples an enlarged ensem- 
ble that also contains the two-point-functions at all possible separations. The "worm algorithms" 
that can be formulated in the loop gas representation of such an ensemble typically show hardly 
any critical slowing down. In addition the non-standard formulation often makes improved esti- 
mators for key observables possible. The basic idea behind worm algorithms is very general and 
can be applied to a large class of statistical models. Efficient algorithms have been developed for 
the one (real or complex) component 4 theory and tested in its Gaussian and Ising or XY-model 
limits as well as at some intermediate values of the coupling [Q, BL H, Bl ^, R]. More complicated 
systems like nonlinear 0(N) — o models [Bj, CP(N — 1) models [P] or two dimensional fermionic 



systems [10, 11] have been successfully simulated as well. Some of the recent development has 
been summarized during this conference [jTJ]. 

In the following we make a first attempt to generalize the idea behind worm algorithms to the 
case of gauge theories. We formulate an algorithm that samples the (generalized) partition function 
of an Abelian U(l) gauge theory. 

2. Strong coupling expansion 

We start from a Wilson plaquette gauge action 

S[U] = -p £ Re[U(x,ii)U(x + fi,v)U- l (x + v,n)U- l (x,v)} , (2.1) 

where x is a site on a D-dimensional hypercubic periodic lattice of extent in direction fj. = l...D 
and U (x,/J.) £ U (1) denotes the gauge field on the link that connects site x with site x + fi. We use 
lattice units throughout. 

Observables in this model are products of link variables 

W=UU(x,Li) M , (2.2) 

x\l 

where we have introduced an integer valued external field j(x,n) = ju.(x) € Z. One is interested in 
their expectation values 

v*>-m- <23) 

Z[j] = J DU e- s[u] U j , (2.4) 

where DU denotes the invariant measure on t/(l) on all links. Gauge and center invariance of the 
action causes expectation values to vanish unless the external field satisfies 

d Pii(x) =0, and (2.5) 
X>M = 0, (2.6) 
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where d* denotes the backward nearest neighbor difference. 
On each plaquette one can expand the exponential 

oo / a \ m+n T jn-m +°° 
m,n=0 \ / n= — oo 

The summation variables n(x;/i, v) = n^ v (x) can be defined to be antisymmetric in ju,v. At this 
stage the original group integrals can be carried out leaving behind constraints. The partition func- 
tion takes the form 

{„} \x,n<v J 
with the abbreviation for the constraints 

S[d*n-j}=l\5 d , nvti{x)Mx) . (2.9) 



The representation eq. ( |2.8| ) has been used as a starting point for Monte-Carlo simulations in [ |13| ] 
and more recently in JT4]]. The authors of [|l4|] find that a variant of their algorithm, that is based 
on an expansion of Z[0] of the Z2 gauge theory only, suffers from critical slowing down with a 
dynamical exponent similar to that of a local algorithm in the standard formulation of the model. 

Inspired by worm-algorithms we pursue a different direction. In spin systems the crucial 
algorithmic idea was to consider an enlarged system that allows for defects related to two-point 
functions. Similarly in our case we consider the enlarged ensemble 

2? = Y,p-\j\Z[j\. (2.10) 

{;} 

The sum is over all possible external fields jfi(x) and, as in [Q], each contribution is weighted 
by some non-negative weight p . J? contains now also graphs with boundaries rather than closed 
surfaces only. The partition function can be written in its final form 

^=i( n (2-ii) 

{„} \x,fl<V J 



Expectation values with respect to this ensemble are formed in the usual way 

1 

T 



»]» = 4?E*M( II 4, vW (y3)]p- 1 [a*n], (2.12) 



{„} \X,fi<V 

and for later convenience and in analogy to ref . [||] we also define an expectation value with respect 
to the vacuum ensemble by 



(( ^ ]))0= «5P*„]» • (2 " 13) 
One can relate observables in the enlarged ensemble to those in the original one. One way to 
measure the expectation value of an arbitrary non- vanishing correlator, given by some jn(x), is to 
find some background plaquette-field kn V (x) that solves the constraint d*k = j. The observable 

0[n] = [I In ^ {x)+k ^ {l3) (2.14) 
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is then an estimator for the desired correlator 

(W) = ((0[n])) o . (2.15) 
In practice k will be zero almost everywhere. For rectangular loops the minimal number of factors 



differing from one necessary in the product eq. ( |2.14| ) equals their area. 



3. Monte- Carlo updates 

We use the freedom in choosing the weight function p _1 to restrict the possible external fields 
or boundaries to the subset containing only one non-intersecting loop with zero winding number 
with respect to all torus directions. Such a loop can be completely characterized by the ordered 
cyclic set of sites through which it passes 

W = {x u x 2 ,...,x P ( W) }, x P{W)+l =xi. (3.1) 

Consecutive sites are nearest neighbors xi+\ = x, + /t;. If we extend the possible unit vectors to 
include negative directions: — ju = —(X, plaquettes can be labeled in several equivalent ways 

?fyv(X) = —n-n V (x+p,) = -n^ v (x + v) =7i_ iU _ v (x + /t + v). (3.2) 



To stochastically sample the sum eq. (2. 1 1 ) we perform a sequence of local updates of which each 
by itself satisfies detailed balance. The updates consist of a proposal W — > W followed by a 
Metropolis accept/reject step. 

3.1 Flips 

If the perimeter P(W) is larger than four we pick a random site x ; on the loop and form 
y = Xi + \ —Xi +Xi-\. If y G W the old configuration is kept, otherwise the proposal is to form W 
by replacing Xj — > y and to adjust the plaquette field accordingly 

(x^ + 1 . (3.3) 

Such a proposal is accepted with probability min(l,qni p ) with 

_ 4-,,_ lf ,(*,)+i(ff) p(W) 

qmp ~ i,.^) P (*»y (3 - 4) 

This update changes the shape of the loop without altering its perimeter. It is displayed in fig. [T[ 

3.2 Shifts 

A site x t and a direction v orthogonal to /},• is chosen randomly. Two sites, y = x, : + fx and 
y' = Xi+\ + v are constructed. A proposal is made in either of two cases 

1. If neither y G W nor y' G W , it is proposed to extend the loop by two sites 

W' = {x l ,...,x i ,y,y',x i+ i,...}, (3.5) 
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Figure 1: Flip (upper part) and shift (lower part) updates. The plaquette field in the shaded area differs by 
one unit with respect to the non-shaded areas. 

together with the necessary change in the plaquette field: n v ^ (x,-) — > n v ^ (x,-) + 1 . In analogy 
to the flips we accept with 



<?shift 



(3-6) 



2. If y = Xj_\ and y 1 = x i+2 the reverse move is proposed, i.e. x, and x i+ \ are proposed to be 
removed from the loop. The acceptance probability is again min{\,q s \^). 

Otherwise no move is made. 

3.3 Non-local updates 

With free boundary conditions the alternation of the two updates proposed so far would be 
sufficient to guarantee ergodicity. On a torus however, there exist contributions to the partition 
function that are not yet sampled. This is fixed by another update: 

A plaquette jc,jU, V / /I is chosen randomly. The subsets of sites which belong to the two dimen- 
sional plane through x spanned by jx and v is denoted by Y. The proposal is to change the plaquette 
field according to 1 



and it is accepted with 



for all x 6 T : 



<7pla 



l liv 



-|-r 4 MV (.r) + l(ft) 



(3.7) 



(3.8) 



The acceptance rate for this proposal is tiny unless the volume is small. 

The last update that we have implemented is proposed whenever W is planar, or alternatively 
whenever it has exactly 4 corners. The proposal is to shift the whole loop by one unit in a direction 
perpendicular to its plane. This update can be built up from the elementary shifts and flips, but it 
turns out to reduce autocorrelation times considerably if it is introduced as a separate update step. 

One iteration of our algorithm consists of order of volume many shift, flip and planar-loop- 
shift updates followed by one plane update. The costs are comparable to those of a Metropolis 
sweep in the standard formulation, i.e. O(volume). 

'The inverse move is realized by swapping ji o v. 
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4. Performance of the algorithm 

Abelian gauge theory is particularly well understood in three dimensions. For the model with 
Villain action, which is believed to lie in the same universality class, the existence of a mass gap 



m and confinement at every value of j8 have been proven [15]. The continuum limit at fixed m is 
believed to describe free massive bosons. A recent numerical work with the Wilson action [16] 
showed that the mass gap is indeed well described by 



mL = 5.23(11)- ^87T 2 j8 exp [-0.25277T 2 /3] . (4.1) 

To test our algorithm we use eq. Q4.1| ) to keep the physical volume constant at mL pa 6 while 
increasing the resolution Lja E {8,16,24,32,40}. As weight function we take an exponential 
that depends on the loop perimeter p[W] = exp(6(P[W] —2)), and tune the "loop tension" d to 
an algorithmically close to optimal value. More precisely, the cost to estimate aL/4xL/4 Wilson 
loop to a given precision is minimalized. We measure expectation values of rectangular Wilson 
loops as well as of plaquette-plaquette correlators. The autocorrelation times and a cost-indicator 
(CPU-time x relative error squared divided by volume) of the average plaquette, the string tension 
and an effective mass from the plaquette-plaquette correlator at separation L/A are monitored. The 
same calculations are repeated with a standard Metropolis algorithm, in which the proposals are to 
change the gauge field on a link e"^ — > e'W+ A 0). In this case the number of sweeps per measurement 
as well as the size of the interval from which A0 is drawn are kept at their optimal values. The 
optimal number of sweeps per measurement is determined only on the 16 3 lattice and kept at the 
same value (i.e. 15) for the others. The size of the interval is determined on each lattice separately. 
Expectation values obtained with both algorithms are consistent with each other and a comparison 
of autocorrelation times and costs is shown in Fig. [| 



5. Conclusions 

We have extended the concept of worm algorithms to Abelian gauge theories. Instead of a loop 
gas we deal with a surface ensemble. The worm is replaced by an open surface and a Wilson loop 
plays the role analogous to the worm's head and tail. All standard observables can be estimated in 
this model. In first numerical tests in three dimensions no significant critical slowing down could 
be observed, but larger correlation lengths will be necessary to make a definite statement. On the 
presently available lattices also a standard Metropolis algorithm performs relatively well. 

Abelian gauge theories with different actions, like the Villain model or Wegner's Z2 model 
can presumably be treated in the same way. Whether matter fields can be incorporated into the 
algorithm or whether an extension to non Abelian models is feasible remains to be investigated. 
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